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ABSTRACT 

The basic approaches to image registration are surveyed* 
Three image models are presented as models of the subpixel 
problem* A variety of approaches to the analysis of subpixel 
analysis are presented using these models* 
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This report summarizes studies conducted on Research 
relative to Automated "ultisensor Image Registration* during the 
period April tgr| - April 10;s? by Or. La v een N • <Rn»l <Pr|n c loal 
Investigator), M r . David Lavlne, t'r, Charles Herman and pr, Eric 
Slu*. This work was supported In Part by NASA Grant NAG^-154 to 
the University of Maryland and In part by l.N.K. Corporation, 
Silver Sorina, mg. The main purpose o^ the*e studies **a* To 
review current approaches to image registration rele v ant to 
Landsat images with emohasis on feature matchino and suhplxel 
accuracy* and to recommend research which should be pursued in 
the area C f registration error and subpiael accuracy analysis and 
estimation. 

Or. M.e. Ramaprivan, information Eatraction Branch, Code 
yJZ* NASA uoddard Space flight Center, Greenbelt, *»0 served as 
technical officer for the NASA grant. The princioal investigator 
is orateful to hr. R«mapriyan tor helpful discussion and for the 
references provided by Or* Rawapriyan on Lanpsat image processing 
procedures • 
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This report deals with, two aspects of image registration 
encountered in the processing of Landsat images: feature matching 
anj sub-pixel accuracy. Many recent research efforts in the area 
o* registration have considered the need to transform images into 
a torsi more suitable for accurate registration. These 
transformations range from simple linear filters to complicated 
procedures based on knowledge cased image segmentation. The 
design of optimal transformations is difficult due to the wide 
range of types of transformations under consideration. 

Certain classes of transf ormat ions yield highly simplified 
images which may be registered either as rectangular arrays of 
pixels or as sets of features. Examoles of such simplified 
images are binary eege and sparse point images. Registration of 
edge and ooint images appears to be robust in comparison with 
many grey-level correlation methods. Edge and point image 
registration are important. In Section ii we briefly review 
basic methoos in this area. We also include a description of 
work we have done with one type of feature matching procedure. 

Section III is a brief review of work on estimation of the 
errors inherent in image matching. This material is particularly 
important, since large errors in a matching scheme preclude the 
possibility of extracting subpixel information from that scheme. 
Included in this section is a discussion of a paper describing 
the selection of control points to minimize matching errors. 

Analytical studies in the area of suo-pixel accuracy have 
peer limited. In Section IV we give an overview of some basic 
ideas in the area. Though the analysis of this section is 
limited, we feel that a more detailed analysis of these model 
situations is feasible. 

Section y contains a description of the research we are 
recommencing on edge-based image matching ana the analysis of 
sucpixel accuracy. Three distinct approaches are describe? » 
ranging from e selection of "test" edge features to a theoretical 
analysis of the maximum accuracy attainable in a mathematical 
mouel o 4 the image. The last section, the bibliography, includes 
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some items not cited in the body of the report* for the salt* of 
como leteness » we have included m»ny articles on image matching 
which have indirect bearing on feature matchino or sub-pixel 
accuracy since they do contain ootentially usefull information on 
filtering and image modeling* 
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The recent development of various registration algorithms 
employing some type of feature matching emphasizes the need for 
more detailed analysis of the capabilities of feature matching. 
In tnis section we describe a feature matchina nethod, called the 
Abstract Vector Matchina Method (AV), which has performed well in 
simulation as well as experimentation with aerial imagery. The 
Landsat registration accuracy requirements reouire high accuracy 
on a large percentage of the images processed, T 0 ensure this* a 
reigstration method must be insensitive to many sources of 
grey-level variation such as crop growth, soil moisture and very 
eocal scene changes. Feature matching methoos offer considerable 
promise in providing this robustness. A variety of studies 
performed c,y l.*'.*. Corporation have led us to believe that the 
A , V • method captures many of the desirable orooerties of other 
feature matchina methods while i ncorporat i no some new approaches 
to tne problem, 

Tne basic A.V. method can be stated briefly. First* find a 
set of points of interest in the image, such as road 
intersections or points of high curvature on curves such as 
rivers. Second, setect certain pairs of those points to be 
joined by vectors. The vectors may be real image edges or 
arbitrary oairings of points, which we call abstract vectors. 
Third* carry out the same procedure on a map or a control image. 
Fourth* ♦or each vector in one image, determine which vector in 
the other imaoe could align with ft on the basis of length and 
possibly other factors such as image intensity on the two sides 
of the vectors. Each oossible matching pair leads to an affine 
transf orma ticn taking one image to the other, such a 
transformation is given by a set of rotation, translation, and 
possible scale factors. The parameters definino a t rans format ion 
may be thought of as a ooint in euclidean soace. Clustering is 
performeo in this space on the set of transformations 
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corresponding to matching vectors* The cluster centers are 
recarded as potential registration transformations* Each such 
transformation is then apolied to the vectors in the sensed inane 
a nc a figure cf matching merit is computed. The transf ormat ion 
with the hioh«*st merit is taken as the correct transformation. 

Two features of this algorithm deserve further comment* 
First* the selection of intersections ana abstract vectors adds 
robustness to the detection process* Matchina procedures based 
upon edges perfo r n poorly if t hr edges become fragmented during 
the uetection process* since the ends of the edges become 
difficult to locate. In spite of this f ranmentat ion , the 
intersection of edge seoment s is often preserved, while matching 
could be performed directly on the set of intersections* the 
matching could be performed faster by imoosing additional 
structure on the ooint set. Selecting oairs of ooints, joining 
them by vectors and matching vectors is considerably taster than 
investigating all possible matchings of ooints. Since 
intersections may ge reliable even when edge endpoints are not v 
we do not require pairs of endpoint* to Correspond to endooints 
of a real edc®. 

The second feature of interest in the ala C rithm is the 
clustering of transformations. Each matcnina of a vector in the 
sensed image with a vector in the reference image gives rise to a 
feasible reoistration t ransf orma tion . Clustering transformations 
corresponds to finding those transformations approximately 
maximizing the number of vectors matched. Clustering enables the 
al a orithm to a^oid examination of many incorrect transformations. 
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An image registration algorithm developed at Rockwell 

International corporation is described in CBelsher et al* 19*91* 
The technique is based on matching high contrast edges* 

experiments have beer performed in both the infrared ana ootical 
portions of the spectrum. Since the oattern matcher was not 
testea on Landsat data* it is difficult to extrapolate the 
numerical results to Landsat registration accuracy* In spite of 
this. the methods of analysis and the experimental results give 
some insight into the possible benefits of subpixel 

inte rpolat ion* 

The feature extraction algorithms were designed to maximize 
the peak to sidelobe ratio, defined as the number of matching 
features at the correct offset divided by the maximum number over 
all incorrect offsets* In order to carry out this maximization 
the probability of the existence of a reference feature given a 
sensed feature was determined and maximized for each feature 
extraction algorithm. The expected number of matching features 
at the correct offset is: 

N Prob(SflR) 

where Prob(SflR) is the prooapilty of a match between a sensed Ij) 
ana corresponding reference C R ) feature at the correct offset, 
and N is the number of features in the reference. At incorrect 
offsets, sensed features are independent of corresponding 
reference features* Thus the expected match at incorrect offsets 
is given by: 

N Prob CS )Prob(R) 

where Prob(S) and Prob(R) are the sensed and reference feature 
densities. Thus the ratio of the expected peak to the exoected 
siaelobe is 

WS> - Prop ' s ) pi£b <ET ■ Pro¥TR7 Pr ° b<R/S) 

where Prob(P/S) is the probability of a reference feature given a 
sensed feature* The reason tor choosina the feature extractor 
maximizing P(R/S) is that Prob(R) can be computed from the 
reference image* 


( 6 ) 


ORIGINAL PAGE 18 
OF POOR QUALITY 


Edge features were obtained by median filtering to remove 
noise spikes* followed by the application of a gradient operator* 
Further processing Is done to remove the effects of local 
variations In grey scale density* and finally a template 
extractor Is used to oDtaln the endpoints* Spot features are 
al so extracted • 

Each Mature matching leads to a possible matching of a 
reference to a sensed Image* For N such possible matching we 
denote the a posteriori probability of offset 1 being correct 
given the sensed image by p(i/S), i=i»*.. t N* using Bayes theorem 
it is seen that 

P<i/S> - P ( ^f) (s/i) 

In general, flight information can be used to estimate p(i>* and 
p(S) Is common to the evaluation of all offsets* so it suffices 
to compute p(S/i)* 

The sensed features* * at the correct offset are assumed 
to take the form 

Sj * Rj ♦ 

where is the reference measurment and H is noise* The 
measured features are mostly lengths* From the above 
decomposition of Sj * the following expressions for p(S/1) 
results: 

plS/i) » p(U> • p(S - R) 

where 7i*"S» and "R are vectors with coefficients aj * Sj and Rj 
resoect ive ly* Finally* assuming the errors are independent the 
following expression for p(i/$> may be derived: 

logCp(i/S>3 * logp(l) ♦ ? logplSj - Rj ) 

Some modifications are required to handle problems with very 
small probabilities* A distributional form for p must be assumed 
with reasonable choices including the Gaussian* laplacian and 
Cauchy densities* 

We may think of p(1/S) as a correlation surface which is a 
function of offsets 1* The surface is approximated by a 


quadratic surface using interpolation on a 3x3 pixel area or a 
>x5 pixel area* Registration was performed on 3-0 imagery using 
fifteen images* The root -me an -square registration error with no 
interpolation was L*3 pixels* Using subpixel interpol at ion the 
average error was reduced to I* 9 pixels* In only one case was 
the error increased using subpixel interpolation* 
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The tr»ue-offs between correlation and feature matching have 
not reen studied adequately* Ratkoyic C1v?93 outlines basic 
issues in assessing this trade-off and describes a hybrid 
correlator combining desirable features of both methods* Feature 
definitions are given and the relevance of feature 

chjrac terist ics to correlation performance are described* 

Accuracy* as described in Ratkovic's work, refers to the 
aicth arouna the true correlation peak in which a correlation 
point is likely to lie, rather than the accuracy attainable 
through interpolation of the correlation function. Three basic 
issues are discussed* Hom ooes one quantify accuracy? What 
factors affect it? what methods improve ft? It is assumed that 
„ two st.ne matching process is available in which the first 
sto?e determines an approximate match, while the second stage 
improves local accuracy of the registration. 

Based on limited experiments, it appears that the three 
primary factors affecting registration accuracy are the number, 
uize, and mean intensity of homogeneous regions in the image* A 
homogeneous region is defied to be a set of soatially connected 
*,ix*is which posses the property of at least first-order 
stationarity and oossible second-order stationarity* Several 
experiments were conducted to determine the sinnificance of 
homogeneous regions* First, a 20 x 20 reference scene was 
generated using a Gaussian distribution with zero mean and 
variance one. The center 10 x 10 inane was designated a sensor 
scene and correlated again s t the reference scene* ?o introduce 
the effect of nonhomogene ity , a bias value was added to all 
values in the left half of the scene* Correlations were 

performed for offsets M = 0,1, 1C, 100. »s the oias level 

increased, the correlation function became prooressively flatter 
near its ce»k* 

further experiments were conducted with CRTS satellite data 
from several tyoes of areas (agricultural, mountain, suburban, 
„nc cesert). The images were divided into homogeneous regions 
unu correlations were performed in two ways. First the ordinary 
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correlation of the sensor and reference imaqe was computed* 
Next, each pixe* In a homogeneous region replaced by the mean 
value of oixel intensity in the region* and cerrelation was 
p»rformeo* Little difference was observed usinq the two tyoes of 
correlations. r rom the above experiments* as well as others* 
rtatktvic infers that the presence of large homoqeneo js reaions of 
*i jnif icantly varyinq mean intensities are the major source of 
peak spreading in the correlation function. 

A "pure" feature matching algorithm was employed using edge 
cetectors to determine region boundaries. Intersection points of 
region oounoar ies were then extracted together and the number of 
lines meeting at such i ntersect i ons was recorded. A matching 
procedure was then aoplieo to the vertices* weighted by the 
number c* intersecting lines to determine the registration 
tr u nsformation. 

A feature matching correlation procedure was developed. The 
images were decomposed into homogeneous regions and each region 
was indeoenoently normalized to zero mean and the resulting 
imdots were correlated* The orimary difference between this 
<*ethod and the "pure" feature matching -sethor* is that in the 
former* region matches are weighted by their si*e. while they are 
not in the latter. 

The final matching algorithm prooosed by R*t kovi c is a 
nytrid between correlation and feature matching. In fhi s 
algorithm the reference image is segmented into homogeneous 
regions. *t each displacement between the reference and the 
sensed image* the sensed image is segmented identically to the 
reference Image* and corresponding regions are normalized and 
worrelated. The correlations in the individual regions are added 
one the registration t ransf ormat ion is computed from this 
correlation function. Y he hybrid correlation algorithm avoids 
the difficulties of reliable feature extraction in the sensor 
image* while retaining the advantage of correctina for 
f'oco.j^nfous regions* 

Feature matching and the hybrid algorithm both essentially 
nigh-oass filter the inane to remove the effects of varying mean 


ORIGINAL PAGE 18 
OF POOR QUALITY 


intensities t»»tween regions. This high-pass filtering* which has 
received considerable attention in reeent /ears* is still not 
•elt understood with resoect to its role in the correlation 
process. 

4 conclusion of the ffatkovic study *as t^at* in the presence 
of ocrely local noise* ordinary correlation perforaed best among 
the aloorithms. For regional noise* such as a bias value added 
to an entire region, resulting from such factors as soil moisture 
or crop aevelopment, the cure feature extraction method performed 
pest. The fc^rid correlator give somewhat less accurate 
reiji st rat ion information than the cure feature extractor* but at 
less cost. 
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? hr estimation of translation parameters in inane 

reci r,t rat ion is a special ease of the image registration w ith an 
oftire distortion. Local registration and false acquisition 
rr obab i l i t ies in this more general case were studied by Mostafayi 
unc Smith They determined the optimal reference image 

sh jo e and sire corresponding to a given rotation and scaling 
cetween tne reference and the sense image* The rotation and 
scaling parameters are determined, prior to registrat ion . using 
knowledge of the sensor position and attitude, in the presence 
of rotation and scaling changes, two reference window oroperties 
rust be oalanced. First the window should be large enough to 
obtain a reliable estimate cf the registration parameters. 
Second^ the window should be small enough to orevent the 
geometric (distortion from degrading registration quality, which 
can be caused by poor alignment of pixels that are far removed 
from the center of rotation. 

^ostafavi and Smith draw further conclusions aoout t he 
structure of ideal reference imanes. For a fixed rotation and 
scaling change, they show there is an ideal reference image size 
and shape which minimized the probability of selecting the wrong 
looe in the correlation function. Ther? is an ideal reference 
image size and shape which minimizes the error in regi st r at ion. 
vjiven that the correct lobe is selected. rhe ideal window size 
in the latter case is smaller than the ideal window size in the 
former case. 9rth the orobability of false acquisition and the 
local reoistration error aepend cn the difference matrix I-A 
where I Denotes tne identity matrix and A defines a rotation and 
scalt cranoe. Before discussing the significance of these 
conclusions* we present the assumptions of the model. 

The reference image. I. (x) * P(x). where x = fx^.x^) is a 
function of two v aridble s *her« the origin is Assumed to be the 
reference point. The sensed images* I. Cx) is given by 

I s (x) • P d (x) + N(x) 

( 12 ) 
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where P<j(x) is a 3 eomet r i cal ly distorted version of the reference 
inane and Mfx) is additive noise* It is assumed that P(x) and 
• if*) are uncorretated* zerp-mean Gaussian random processes which 
have soatially s hi f t'invar iant » second-prder statistics* In 
addition the cross-correlation is assumed to t>e Gaussian* The 
patterns arp assumed to he differentiable in the mean-square 
sense. The geometric distortion is given by 

■ P(A X + c o> 

where A represents a rotation ano a scale change and tp is a 
translation. The basic measure of correlation acquisition 
accuracy is the expected value of the correlation process of the 
true translation divided Dy the standard deviation of the 
correlation process far from this point. The local registration 
accuracy is a function of the gradient of the correlation 
function at t^ • For a correct registration. this gradient should 
be zero. 

A basic difficulty in applying this work lies in the 
assumptions made in the analysis. The expected value of the 
correlation function 3t t^ » and the standard deviation of the 
correlation *ar from tQ are computed independently* Thus a more 
realistic but less tractable problem is to look at the value of 
C(t Q > and ^ = C max C ( t ) | t far from tg > for each realization 
cf the Gaussian picture process and determine the probability of 
the event > C(t 0 >. Similarly, the value of the gradient of 

Cft) as a measure of reoistration accuracy is difficult to relate 
directly to the expected registration error in Pixels* 
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The variance of the translation error In registration Is one 
measure of registration performance* Two procedures for 
estimating this variance CMcGillem and $vedlow 197 J will be 
described in this section, while the variances compared using 
these methods may oe less than a pixel. the authors assume that 
the log likelihood function for the image to image displacement 
can be modelled by a second order polynomial. This assumption 
coes not appear to be derivable from more fundamental properties 
of their model* nor i s any experimental evidence for its 
plausibility provided, As a result of this assumption, the role 
of various correlation function interpolation schemes in 
registration accuracy is bypassed. thus limiting the 
effectiveness of these models for error variance. Though the 
relevance of these variance estimation schemes is currently 
limited by the quadratic model, the basic iaeas of these methods 
will be presented due t° the possibility of future improvements 
on the methoo. 

jn both methods, it is assumed a model for the underlying 
image is known. Method 1) makes the following additional 
as sumptions z 


1) The sensed image consists of additive noise 
which is independent of the model. 

2) The joint probability density function of the 
noise is Gaussian. 

3) The a priori distribution of the translation 
parameters is uniform over the range of 
interes t . 

4) The variance in the x and y directions may 
be modelled separately. 

5) The final result is dependent upon a targe 
signal to noise ratio. 

The Likelihood function of the displacement parameters is given 
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where 


A(T r = likelihood function of the 

x y 

the displacement parameters T and T 

X J 

P m <T^ » T y ) = density function of the parameters* 

T and T * ci v en the known signal 

* y 

P m (T^»Ty)(f) = conditional density of f(x,y) 

given m( x*y *T_ ,T„) is present 

A y 

Pq f f ) * conditional density of the f(x f y) given 
m (x t y f T x ,T^) is absent 

m(x*y*T » T ) = known signal as a function of 

x y 

the spatial coordinates and the 
displacement parameters 
F(x t y) = a (x t y) * n(x«y) * received signal 
n(x t y) * additive noise 

The variance of the error in the x direction translation 
■ ill now be estimated* Denoting the estimates of T x and T y by T x 
and Ty* first expano the log-likelihood in a Taylor series about 

T . and then u se the fact that at a maximum of the likelihood 

x 

function* the first order derivatives are zero* The following 
expression is obtained: 

In A(T T)~?ln A(T T ) + 1/2 ^ A(T *» V (T _ $ > 

7 7 5T X 2 * 7 

The variance in the x direction is then given by 

r.2 ,A A „ 

2 9 In A(l' x , T y ) _i 

^ ■- [ — zr* — 


cx pand in. 


A (T , T ) it :an be shown that 

x y 


J*. - ZZ 0 , W V 3m h (T x’ V 

V gh — — 
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Here, the are terns in the inverse of the covariance matrix 
uf the noise* g and h represent index pairs gsCx^yj)* hsfx^y^)* 
and the sun is taken over all indices* The above expression can 
be simplifiec to yield a variance estimate in terns of the 
effective bandwidth and signal to noise ratios in the x and y 
di rect ions • 

The second method assunes the noise is additive and 
independent of the signal; the noise spectrum is known and a 
matched filter is used* The basic model assumes the received 
signal is a sum of t«o output signals; the model convolved with a 
filter and the noise convolved with the same filter* First the 
functions representing the filtered model and the composite 
filtered output signal are expanded in Taylor series about the 
true registration position* It is assumed that the filtered 
outputs are maximized at the true registration position* 
Applying the corresponding constraints on the derivatives* an 
expression for the error variance is obtained in terms of the 
filtered output noise and the filtered model output* By 
selecting a matched filter* a variance estimate similar to that 
of methoo one is obtained* 
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T he spatial distribution of ground control points (GCP) can 
affect the accuracy of image registration. Clearly the number 
ano location cf SCP's is constrained by the scene under 
consideration. but subject to this *.onst ra int , a registration 
proaram must still select a set of ground control points. This 
selection orocedyre s h ou^ d Choose Control points so as to 
minimize the reoistration error. in a study of thi s Dro ble w » 
Orti Cf9ri7 developed an expression for the mean sauare Landsat 
MSS registration error as a function of GCP location s# A 
minimization procedure is then apolied to determine the ootimal 
GCP ci st r ibution. 

Ortizs procedure estimates the attituae and altitude of the 
satellite *rom a set of SCP's. All other mapoing functions 
relating the raw and corrected images are assumed to be known. In 
adoption it is assumed that within a Landsat frame the three 
attitude ancles can be adeauateiy described by cubic polynomials 
o f time ana the altitude by a linear function. Thus there are 
fourteen coefficients to ue estimated. The error criterion is 
the averaae error in these coefficients. 

The geometry of the transformation problem will now be 
describee. All deviations are assumed to be small. The attitude 
a n o altitude deviations from the nominal values are given by 

x ■ <J> + ktanG(p) 
y ■ w[l + tan^GCp)] + htanG(p) 

where p is the CC° , 's oixel number within the scan line measured 
from its center. &(o) is the Corresponding scam angle, zero at 
thr center of the scan line; <p ,n , and x are the attitude 

uncles, citch, roll, and yaw and n is the relative altitude 
deviation, all measured at the time when the GC D was seen by the 
satellite* Here x and y reoresent the GCP control differences in 
position along and across the orbital track on a plane tangent to 
the earth through the nadir point of the imaoe. The x and y 
coordinates represent differences in positions between the GC P in 
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map coordinates and the GCP position obtained from transforming 
the raw image coordinates without consideration of attitude or 
altitude deviations from tne nominal values. 

The attitude angles are then aoproxiaated by third degree 
polynomials using orthonormal Leoendre polynomials and the 
altitude is apo rox i mat ed by a first degree polynomial* The 
expressions for x and y together with the Legendre oolynomial 
expansions yield the following matrix equations for deviations 


lx »y >; Ci=1*.**»n) corresponaing to GCP's c and c 

i l x« 


x - W x • c x + C x 


Y - W.. 


C„ + t 

y y 


and F are vectors reoresentino the errors associated with the 

a y 

control points. w and W are matrix functions of F and the 

x y 

Legendre polynomials entering in the expansions. From the above 
matrix equations* least-mean-square estimation of c and c can 

* y 

te computed* These estimates can then be used to compute the 
variance of the oropagated registration error. After a somewhat 
involved derivation of the registration error as a function of 
GCF location it is determined that GCP"s should lie near the four 
corners of the image and in some areas o* the left and right 
edges of the image where the oositions «.re a function of the 
number of control points* 
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Three approaches to the estimation of edge locations to 
sub-pixel accuracy are presented in a recent paper of Hyde and 
Davis C19823. In two of the approaches* the image is assumed to 
be a "grey level” image in which the intensity at each image 
point can take on a range of values* A set of pixels is assumed 
to have been identified as belonging to an edge which is assumed 
to separate two regions* The regions are assumed to have either 
constant* but Different* grey levels* or grey levels given by two 
different but known Distributions* 

for each pixel ana each possible line passing through the 
pixel* the area on each side of the line is computed. The grey 
level of an edge pixel is assumed to be a weighted linear 
combination of the random variables giving the grey levels in the 
two regions* The weighting is by the relative area of the pixel 
belonging to each region* The edge pixels' grey levels are 
assumed to be inoependert* If we parameterize the edges by polar 
coordinate ( p» 0) * then the probability density function for the 
observed grey level in an edge pixel can be calculated* This 
density is parameterized by (p*0)» since we get a different grey 
level mixture Depending on how the edge intersects the pixel* 

The estimation of the edge parameters* ( p* 0 ) * can now be set 

up as a maximum likelihood or least squares procedure* This is 

straightforward to set up since the intensities are independent 
random variables with known densities* The likelihood estimate 
is difficult to compute numerically and no tests were performed 
using this approach, 

A simpler approach to the problem was also investigated* 
The edge was estimated to be the best least squares fit to the 

centers of the edge pixels* No use was made of grey levels in 

this approach, and no optimality was claimed* A refinement of 
all the above methods was also studied. In this refinement, 
edges were fit to each of three consecutive points on a line, 
rather than globally fitting as described abo v e* 
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The Methods Mere tested on ten Landsat inages in which the 
edges consisted of agricultural field boundaries. The edge 

pixels were selected using ground truth rather than an edge 
detector, as would be neces s ar y in practice. Bands 1 end *• were 
selected, since they tend to be different and independent. 
Results are reported for least squares fitting using grey levels 
and least squares fitting using pixel centers* Both local (three 
point) and global fitting experiments were performed* For global 
fits, ground truth and estiaation results for p and were given* 
For local fits the rat errors in 9 and p were reported* The rat 
errors in p were less than one in all cases and the ros errors in 
@(in degrees) ranged froa about *3 to 12, with results clustering 
about 1*5* The fitting to pixel centers M as a$ accurate as the 
edge fitting using grey levels. Thus the use of grey levels did 
not appear advisable, due to the additional coaputational 
coaplexity. 

The experiaental results of this work are difficult to 
assess because the accuracy is given in terns of p and 0 rather 
than in terns of pixels. 6lobal Measures of registration 
•ccuracy are not readily conputable iron the given data. This is 
not a criticise of their Methods, but rather indicates the need 
for further work before their study can be adequately assessed. 
No analysis of the probabilistic properties of the estiaated 
registration transformations are given in this work. We believe 
that an analysis of this type is feasible for Methods similar to 
those used in this study. Such an analysis could provide a 
useful basis for the evaluation of algorithas for registration 
with subpixel accuracy* 
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T^o basic approaches to sub-pi*el registration accuracy are 
discussed in a recent survey Cfolfe and Juday 1°8?3 on 
registration. In the first approach, the correlation between two 
images is interpolated and its maximum is taken to he the offset. 
In the second approach, patches in the two images are matched to 
the nearest pixel using correlation, and a mapping function is 
fit to these matches. The success of this method is based on the 
assumption that the errors ir. local patch fits are random and 
tend to cancel. 

Various approaches have been taken to the interpolat ion of 
the correlation function. Fitting of lower order polynomials to 
the correlation function is one common method. Both bivariate 
and fourth-order polynomials have been used for this purpose. 
Fitting is commonly done over a small neighborhood, up to 5*5. 
Ihe order of the p.>. ,r.omial can be allowed to vary. The quality 
of the peak in an interpolated correlation function can be 
evaluated in terms of the curvature at the peak and the heioht at 
the peak. 

Several variations to the above approaches have been tried, 
jivariate gaussians rather than polynomials have been used for 
fitting; elliptical cones have also been studied. The centroid 
of the correlation surface has also been considered as an 
estimate of image offset. As a final approach, the image offset 
can be estimated by the phase function in the correlation 
transform. The offset can be directly calculated from the phase 
tf the correlation function is symmetric, which unf ortunate ly may 
not oe the case. 

The above methods are not based on any theoretical model of 
the correlation process. As a result* t**e accuracy of 
interpolation depends largely on the surface fitted to the 
correlation. r xperinental studies do not provide an adequate 
basis for a clear comparison of these methods at this time. 
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Interest In M^orafng loage registration using feature 
•atchlng Methods Is grossing In the loage processing conounity* 
The teaporal stability of edge locations Is considerably better 
than that of grey levels for a vide variety of scenes. 
Correlation of Inages filtered for edge enhaneeeent May be 
regarded as an Internedlate step betveen traditional correlation 
Methods and More syebollc feature Matching* In this section ve 
exanine somm issues in edge based loage Matching* Ue think uork 
along these directions could shed considerable light on the 
nature of sub-pixel accuracy and the desi9n and evaluation of 
algorithns to achieve it* 

The study of sub-pixel accuracy in registration nay 
profitably be divided Into tvo areas of research* first* how 
ouch inforoation about sub-pixel accuracy is contained in a pair 
of inages* Second* how algorithns can be designed to achieve 
sub-pixel accuracy* In part A of this section ve deal with the 
first topic* and in part B with the second* 


» 
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The Information content of ed*e pixels Is of considerable 
Importance for several reasons. Knowledge of the maximum 
registration accuracy attainable gives a benchmark by which to 
evaluate algorithms, It can reveal situations In which It Is 
fruitless to pursue the sua*p1xel problem, and It may give rise 
to tore accurate registration algorithms. 

The design of models for feature based matching Is at a much 
more primitive level than Is the modeling of grey scale Images. 
A variety of random field models are currently In use In the 
study of grey scale Images. Though there Is considerable 
question as to the usefulness of these models, they have proved 
useful *s * study point In the analysis of the relgstratlon 
pro ces s . 

We confine our modeling discussion to binary edge Images. 
Various arguments can be made for and against such a restriction. 
On the positive side, such Images are easily attainable through a 
variety of edge detection techniques. The simplicity cf the 
Imane greatly reduces the information content of the picture, 
thus allowing more extensive processing to match Images. This 
simplicity also offers hope that the situation can be analyzed. 
On the negative sioe, the relationshio between th* binary image 
and the original Image Is no loneer clear. Unless reliable 
information on the accuracy of edge detectors is available, it is 
difficult to evaluate the accuracy cf the binary Image 
reg 1st rati on. 

The objection raised above to the use of binary images for 
reg istration can be softened by several factors present in the 
Landsat context. Registration procedures currently used with 
Landsat images can register to within a pi«el on a large numter 
of scenes. Current attempts at suo*pixel accuracy separate the 
registration procedure into rough and fine registration, where 
rou:h registration is cone ty correlation and fine registration 
is done by interpolation. i similar apnroach can be aoopted with 
geometric matching. 

The exoense of reliable edge detection can te reduced 
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considerably if prior knowledge of apcrox iaat e edge location is 
available <. To use such knowledge* a rough registration is 
recessary* One approach to regie ration along these lines is to 
perform a rough registration usine correlation* ana an eege 
database to perform edge location in a small number of areas 
where edges are known to be present* Knowledge of edge 
properties, such as spectral si a nature* together with the edge 
location information available from the rough registration* could 
greatly enhance edge accuracy* Edge points such as roads could 
also be sought* in this constrained edge searching environment* 
estimation of edge location accuracy may be plausible* 

The prior edge knowledge allows us to restrict ourselves to 
a very small number of reliable edges* In this section we 
restrict ourselves to a single edge or a pair of patallel 
straight edges* This restriction is imposed because it is the 
natural first step in the analysis of binary edge images* In 
practice* this condition is enforceable by throwing away all 
other edges. 

we now come to the fundamental auest ion of this section* 
Given a set of pixels representinc the digitization of an eege* 
to what extent is the exact location of the edge cetermineo? 
This question requires some amplification before it can be 
formulated precisely* Do we view a real world edge as having 
thickness cr as an ideal line segment? which pixels represent an 
edge? If an Ideal line segment Intersects a pixel very near a 
corner* should that pixel indicate an edge point? Should false 
edge points be allowed? Though these considerations are somewhat 
obvious* the specific model adopted can greatly affect the course 
of analysis* 

we first consider the noise free model* In this case* each 
pixel whic* is set to one must contain at least part of an edge* 
All ether pixels are set to zero* 'ote that we have not said 
that each pixel which contains part of an edge must be set to 
or.e, only its converse* Tnus we allow 4 or the possibility that 
our sampling mechanism requires that an edge be at least a 
certain distance into a pixel before the pixel is set to one. 
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As a first step toward a more careful formulation of our 
basic aoestion* we ask what it means for two line segments to be 
close* There are many definitions which could be used* but we 
must recall that whatever definition is adopted must be useful 
for re 3 ist r it ion* viewing tine segments <i s sets* the standard 
definition of distance oetween sets (minimum distance between 
pairs of points* one ir. each set)* appears useless* oecause cf 
tne possible wide variation in line position among lines close in 
this sense* Similarly* viewing lines as pa ramet eri jed by polar 
or cartesian coordinates* we can define closeness in terms of 
nearness cf the corresponding parameters* This also *an result 
in considerable separation between limes unless line segment 
lengths are -'tplicitly considered* 

Several feasiole definitions of line segment clo$eness can 
be given* we now describe one appealing measure* Let and 
denote two continuous (not digital) line segments* For each 
pixel xg L^ » let d^(x) cenote the Euclidean distance from x to the 
nearest point in L2» Similarly, for ye L 2 » lot 32**' denote the 
distanc* from y to »he nearest coint in L^ . Let 

^ djlx) * 

x £L^ 

D 2 = rcax d £ <y ) » 

y £ 4 

ana C = m a x ( D ^ » D 2 ) • 

Then 0 is a measure of line closeness whicn guarantees that lines 
never get tar aoart in the sense tnac no point on either line is 
very far from the ether line. For a pair of line segments giv^n 
in cartesian coordinates, formulae for this measure can be easily 
de r ivea • 

Several possible definitions of closeness have been 
presented. Of these the last was ch a ra c t er i z ed as feasible for 
re$ ist rat i on. This statement reouires some justification. 
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„Uj.Dc$e we hav» ooserved a set of oixels reDr*s«ntin3 on edge, 
assume furthermore 1 1. . t any two real edges* which could have 
caused this set o* pixels to he set to one, can be oroven to be 
clcst in the aOove sense. What can we say about the true 
location c* the line in the image? Recall that the segment we 
ore observing digitally is known to correspond to a straight 
segment * e.n., a section of road between two intersections. 

We wish to determine a rigid mapping from the road segment 
into the image. Suppose we require the maoping to send the road 
segment center into the center of one of the possible lines going 
through the observed pixels. Suppose we further require the 
mapping to send a whole interval surroundino the road center to 
an interval surrounding the line seament center. Then it can be 
shown that if any two feasible lines are sufficiently close by 
cur cistonce *ess u re h * the n the maximum registration error can 
ce mode less than a pixel and in fact the error can be bounded. 

Cven in the simple situation described above* the hard oart 
of the task remains. Given a set of pixels how does one 
c«te.mine the maximum distance between feasible lines. A related 
question which is of vital importance in this area is the 
following. 1$ it tru 0 th a t for most tine segments which could be 
fluCtd on a sampling grid, the resulting $et o f oixeis is rigid? 
ty rigid, we mean that the set of feasible lines f or that segment 
„ r e dost. A further question of greater practical importance is 
.hut lennth and directional prooerites of a line segment 
guarantee rigidity? Finally, how can a set of line segments be 
selected to guarantee rigidity if single segments are unreliable? 

The answers to tre above questions are not yet available, 
put we are currently investigating these questions under another 
research orant TnaSa Subcontract L2C0C813, *e briefly describe a 
method of attack which we consider to be fruitful. A line is 
ri^i,. if slight changes in its location will result in a chanoe 
in the set of pixels taki.ig the value one. We are currently 
attemetino to show that certain relations between the line length 
anu the line angle guarantee rioidity. Rouohly speakino, the 
conditions guarantee that the line runs near the Corners pf 
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several pixels. Furthermore, the oattern is approximately 
alternating in that the line goes near a lower corner of a pixel, 
then near ,n upoer corner* then near a lower corner, etc. These 
corners may separated by several pixels. 


This concition can oe formulated in terms of certain basic 
functions occuring in elementary number theory. The statement 
that th« condition holds most of the time is related to some 
basic results in probabilistic number theory. 
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The sub-piin edge location information ©resent in a binary 
edce has been disucssed. Ut no,, consider the Question of how to 
extract this information* ’he problem may be subdivided into two 
narts* First* for a given binary edge* how can the set of 
feasible candidates for the true edoe be determined* Second, hov 
should one candidate be selected from this set as an ootiaal 
estimate. The detailed study of these questions is currently 
underway as part of a NASA grant CNASA Subcontract L20QQ813. 

The determ ination of the set of feasible edge locations 
based on an observed binary edge is a complex geometric problem, 
even in the optimal case where all observed edge pixels really 
belong to an edge* In this case, one could rouohly fit an edge 
10 the pixels using scm e method such as least squares* and then 
^erturb the edge till it oasses through the desired pixels* If 
H ixel errors are allowed, then the edge might be constrained to 
go through at least a certain percentage of the pixels* 

Cnee a single feasible line is located, the line can be 
perturbed to estimate bounas on the edge location. The precise 
manner in which these perturbations should be done must be 
determined by a detailed geometric analysis of the oroblem, which 
is currently under way* If reliable bounds can be found on the 
variation in edge location, then several aoproaches to the 
selection of an ootimal edge can be adopted. first, one can 
define a distribution on edge locations, and compute the expected 
edge location. Due to the lack of information as to a reasonable 
distribution, a second aoproach may be preferable. in this 
approach, the "most central" line is located* This line is that 
feasible line which has the smallest maximum distance to the 
ether feasible lines* This can he made precise in several ways* 
out the comout a t iona l crocedure for locating this line will 
depend heaoily on the geometry of the set of feasible lines. 

The feasibility of these geometric methods cannot be 
assessed at this time. They offer oromise because they are based 
on a model of the relationship between binary images and edge 
locations. The major experimental study necessary to tune the 
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mouel is an analysis of the relationship between observed edge 
pixels and real edges* The primary quantity to be estimated is 
the percentage of detecteu edqe oixels which are true ecge 
eixrls. ’his quantity mioht be estimated as a function of eoge 
snrte* length and content (road, field, etc*)* 

The study of edge location coulu be extended to more complex 
figures such as intersections, or strips such as wide roads* 


0RK3M- & 

OF POOR QUALITY 


( 29 ) 


ORIGINAL PAGE IS 
OF POOR QUALITY 




v.i icausis 2i iiaacJifias 

a . s^sLiisycii 

In analyzing the problem of attaining sut-oixel accuracy* we 
■»ust first ask whether, under general conditions* sub-pixel 
accuracy is acMevacte with a reasonaole aegree of confidence; a 
negative answer renders all other Questions moot, Assuming a 
positive answer, the analysis can be divided into two parts: 3 
stu "y of existing algorithms, and the development of new 
algorithms which satisfy certain conditions of optimality. 

In performing the above analysis* one of three (not 

necessarily disjoint) methods can be applied. The first is a 
t r i a I- a nd- e r ro r -t»y - exp e r i m ant a t ion approach* in which "real*' data 
are used *or the experiments. This method is the least 

satisfying, since without a deeper un de r s tand ing of the problem 
(which is exemclified tv a model) real analysis is hardly 

possible. However, it is reasonable to test algorithms on such 
real data. 

In t h " second m c thot, a .nccel is h y po t n es i z eu • The 

consequences cf the mocel are derived v i J simulations* An 

important advantace of this method is that relatively complicated 
moolis can be studied. 

The tnird net hoo also involves a model* however this time 
all results are cerivec analytically. "o guarantee t r act ab i l i t y * 
only simple mooels can be considered. Sometimes the model must 
be simplified to the extent that the very dualities of the 
pronlem which lead tc fundamental difficulties are assumed out cf 
existence. 3 ut if one feels that the model does capture the 
essential f°atures of the problem, then this analytic metheo is 
superior tn the others. It can lead to deecer insights and 
superior algorithms. 

The work currently being cone is restricted to the last two 
methods, with an emphasis on the last one. The model used, and 
the analysis performed, are extensions rf the work done oy Novak 
[19-13, h e considered the problem of finding an edge, of xnown 
shape, in a gr e y level picture. Tnis is essentially the same as 
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•atchina a binary image to a grey Level image* Novak assumed the 
pixel grey levels are statistically independent* and within the 
homogeneous regions separated by the edge* are identically 
di s tri buted • 

To determine the location of the edge within the pixel* 
Novak assumes that the edge pixels have the same distribution as 
the other pixels* but with a mean which is a convex combination 
of the means of the pixels in the two homogeneous regions* The 
weights reflect the percentage of the edge pixel residing within 
each homogeneous region. Using this model Novak finds the mlE of 
the weights* which determine the location of the edge* 

Because this method requires a knowledge of the means within 
each homogeneous region - these means are unknown - Novak 
recommends fitting a biquadratic surface to the nine correlation 
points centerea at the point with maximum correlation* The peak 
of this surface determines the sub-pixel match point* 


3* E§£22E£E2£S ifirk 

Using computer simulations we recommend determining whether 
sub -pixel accuracy is attainable* Because this problem would be 
studied via simulations* rather than analytically* a complicated 
underlying moael can be assumed* Rather than two homogeneous 
regions* we coulo allow many such reqions* with grey levels 
varying linearly over each region* Edge pixels would be assumed 
to be convex combinations of the grey levels of the neighboring 
pixels, rather than simply their means* We would also allow 
statistical dependence of pixel grey values within regions* 

maintaining inaependence from region to region* Note that we are 

allowing only translation errors between the two images* 

Assuming favorable results from the simulation experiments* 
we would proceed to analytically study current matching 

algorithms* The model used is similar to the one described 

above* the difference being that the covariance between pixels of 
a region will take a simpler form in the present case* 

The algorithms would be analyzed asymptotically for four 
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reasons* First, asymptotics are easier to deal with, hence more 
comolicated mooels can be treated* Second, images contain enough 
pixels to make asymptotics reasonably valid • although how many 
are "enough" should be looxed into* Third, asymptotics are often 
independent of underlying distributions (e*g*t the Central Limit 
Theorem), hence results are valid for a large class of 
distributions* And fourth, it is an axiom of statistics that 
estimates should behave properly in the limit, hence the results 
of matching algorithms should approach the correct point as the 
size of the image increases* 

Lastly we would seek to devise new algorithms based on the 
insights gained in the above work* It is clear that in general 
biquadratic interpolation does not behave properly in the limit - 
in fact this remark probably applies to any polynomial surface. 
The appropriate interpolating surface should be determined by the 
aut ocor rela t ion structure of the image, and in fact should peak 
at the same point as the likelihood surface over the image* It 
should be kept in mind that implicitly we are assuming that we 
can only see the (continuous) correlation surface at a discrete 
number of points, and that a second order Taylor series expanded 
about the maximum point is valid over the entire 3x3 neighborhood 
of that point. 

This assumption is not always valid, and in fact there has 
been some confusion in the field about the relationship among 
least squares estimates, maximum likelihood estimates, and 
correlation when the underlying distribution is Gaussian, It is 
assumed they are the same. This is not true in general, and we 
recommend a thorough study of their relationships under various 
con dit ions* 
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We noy oescrioe a general stochastic model underlying and 
extending the framework of «ostafavi and Smith T7gl and allowing 
mathematical formulation of one fora of the problem of sub-pixel 
accuracy of image registration* We suppose that over a region of 
interest there exists a reference grey level field Z R <_x> Ot is 
a point <x»y) in the clan?), including noise* which may be 
modelled as a (strictly) stationary random field* i*e«* a 
collection of random variables all of whose joint distributions 
remain invariant under translation of the olane origin* The 
remotely sensed grey lev*l field is modelled as the sum Zg (x) = 

Z D *x> ♦ Mix)* where 2^ <x) is a distorted form of ^ (x^ ) 

(wostafavi and smith considered the case of affine "geom^tr ic** 
cistortion Z D (x) * z R (A_x * where A is a 2*2 matrix) and N <x) 
is a stationary noise field assumed independent of (Z D (x^) *Z R (jt ) )• 
•e further assume that Z and Z are jointly stationary (have 

D K 

jointly translation-invariant statistics) and that the fields Z^* 
Z R * M have means Z (i.e,* means have been subtracted away) and 
known (or estimable) covariance functions 

fL (t) * Co vfZ_ ( x ) , Z_ (x ♦ t)> (for x»t in the plane) 

0— R — R — — 

P (t) = CovfZ^lx) *Z Q (x^ ♦ £) > 

° 01 lt) = Cov« R (x) f ip * t_)> 

P 2 (_t) = Cov{M(j|) ,N(jc ♦ _t)> 

(which do not deoend on x because of stationarity) • 

Moreover, : (•) is (due to the digital transmission of 

s 

senscr images) for practical purposes only observed at lattice 
points l Xq t jh, y Q + kh) in the Plane, where (x Qf y 0 ) an unknown 
offset from the origin of the reference image* h is the width of 
one pixel, and j and k are inteaers* In this setting the problem 
of image renistratiom i s to recover (estimate) **o ,y O* * rom the 
observations of 7 R (jc > over a large reference area and of ^ <j» ) 

over a much smaller area * in the plane* The statistics most 
often used for this estimation are the offset empirical 
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2 Z s (lh,mh)Z R ((l + j)h,(m +k)h) 
(l,m)eMg points in window Mg 


where Pg is a subset of (flostafavl and Smith assume all fields 

Gaussian* which we do not* and ao not deal with the discrete 
aspects of pixels)* 

There is extensive probability literature (see references) 

showing that for very large windows W (i*e** in the limit of 

/\ 

infinite window size)* sum-st at ist i cs of the form C.considered as 

j 

processes on the integer lattice in the plane are asymptotically 
Gaussian (as processes* and not simoly as random variables) unoer 
various assumptions on the rate of decay of dependence among grey 
levels in groups of pixels which are widely separated* For this 
reason* our moael problem is to recover ( Xq *yQ ) as best we can 
from knowledge of *o ,R Of*l » R 2 an< * observation of Q. j^r all j*k in 
(some subset of ) Kg. Note that we seek not simply to- detect the 
"best" (j*k)-offset approximating ( x^ * y^) b, but to estimate (Xg 
* y ) itself. Cur criterion of performance will therefore be the 

0 A A 

probability that an estimated offset (X*Y) differs from the true 


(x. 


» y g ) by more than a distance (expressed for example as a 


fraction of a pixel wioth)* 
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From the survey of published work on the analysis of 
registration accuracy* presented In the preceedlng pages* It may 
be seen that little has been written on the analysis of subpixel 
registration accuracy* We recommend research be performed on the 
following aspects of subpixel accuracy: 


• A study of the feasible registration accuracy 
attainable using binary edges for matching* 

• k study of subpixel accuracy In Images consisting 
of regions with linearly varying mean grey levels 
using correlation methods* 

• An extension of the random field model of Xostafavl 
and Smith to the non-Gaussian subpixel case* 
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